#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jun  7 15:29:11 2023

@author: jcfq2
"""

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
import shapely.geometry as sgeom
from cartopy.feature import ShapelyFeature
from scipy.io import readsav
import cartopy.crs as ccrs


# In[0]
# set up values


maxint=1.0
minint=0.2
mincount=0
maxcount=3600*100.#3.*3600#hours - many hours due to bleeding of individual pixels across multiple bins
mindate=10
maxdate=30


minrange=0.7
maxrange=1.9
minrange=0.0003
maxrange=0.0018;0.0001
xxrange=6
startx=1

polelat=75
polelong=180




# In[1]
# load data


results='/Users/jcfq2/OneDrive - Northumbria University - Production Azure AD/results/'
planetary='/Users/jcfq2/OneDrive - Northumbria University - Production Azure AD/python/planetary_values/'


idl_file = readsav(results+'jupimage_phase_total_crazytown.idl')
# idl_file = readsav(results+'jupimage_phase_total_crazytown_avgclicks.idl')
# idl_file = readsav(results+'jupimage_phase_total_crazytown_avgclicks.idl')
a=idl_file['a']
a2=idl_file['a2']
b=idl_file['b']
b2=idl_file['b2']

aa=a
aa2=a2
bb=b
bb2=b2


ab=a/b
ab=np.nan_to_num(ab)

   # In[3] set up vogt


# for cml in range(50,329,10):

               


# In[4] 

totalmap=np.sum(ab,axis=0)
# plt.imshow(totalmap)
# plt.show()

totalmapN=totalmap
totalmapS=totalmap
# totalmapN[0:70,:]=0
# totalmapS[90:,:]=0

totalmap=totalmap/np.mean(totalmap)
# totalmapN=totalmapN/np.mean(totalmapN)
# totalmapS=totalmapS/np.mean(totalmapS)

lon = np.linspace(360, 0, 360)
lat = np.linspace(-90, 90, 181)
lon2d, lat2d = np.meshgrid(lon, lat)


crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))


# fig = plt.figure(figsize=(10, 2),dpi=300)
# for i in range(startx,xxrange):

    
#     ax1 = plt.subplot(1, 5, 6-i,projection=ccrs.Orthographic(360-polelong, polelat))
#     mapfull=ab

#     mapmap=np.sum(mapfull[75+i*30:75+(i+1)*30,:,:],axis=0)

#     # mapint=np.sum(a[75+i*30:75+(i+1)*30,:,:],axis=0)
#     # mapcount=np.sum(b[75+i*30:75+(i+1)*30,:,:],axis=0)

#     # mapmap=mapint/mapcount
#     # mapmap=np.nan_to_num(mapmap)


#     mapmapT=totalmap
#     mapmap=mapmap/np.mean(mapmapT)

#     countmap=np.sum(b[75+i*30:75+(i+1)*30,:,:],axis=0)
#     totalmap3=np.sum(a[75+i*30:75+(i+1)*30,:,:],axis=0)
    
#     print(mapfull.shape,mapmap.shape,mapfull[75+i*30:75+(i+1)*30,:,:].shape,lon2d.shape)
    
#     gl = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='black', alpha=0.5, linestyle='dotted', draw_labels=True)
#     gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
#     gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)

#     ax1.contourf(lon2d, lat2d, mapmap,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot') 
    
#     ax1.set_extent([0, 360, 50, 90], crs=ccrs.PlateCarree())
    
#     # plt.scatter(sector2_lon,sector2_lat,transform=ccrs.PlateCarree(),c=sector2_loc,vmin=0,vmax=24,cmap='hsv',marker='.')

#     # ax1.set_global()
#     # 


    # for i in range(1,6):
    #     print(i,75+i*30,75+(i+1)*30)
    
    # for i in range(1,7):
    #     print(i,80+i*25,80+(i+1)*25)
    
    # for i in range(1,8):
    #     print(i,90+i*20,90+(i+1)*20)







import sys
import os
from astropy.modeling import models, fitting

sys.path.insert(0, '../../planetary_values')
from grille_vogt_localtime_dawndusk_all import plot_vogt_slices
from grille_vogt_localtime_dawndusk_all import load_vogt_slices
from grille_vogt_localtime_dawndusk_all import vogt_ofl
from grille_vogt_localtime_dawndusk_all import fit_fieldline_fixed
from grille_vogt_localtime_dawndusk_all import fit_fieldline_scattered
from planet_magfield import Bmag

from shift_magnetic import to_magnetic
from shift_magnetic import from_magnetic


sys.path.insert(0, '../../tss')
from basic.polygonmask import maskeqvalue
from basic.polygonmask import masksubarray



# In[4]
sys.path.insert(0, '../../planetary_values')
   
ofl_lat,ofl_lon=vogt_ofl()
    



# In[5]
brightness=np.zeros([xxrange,36])
brightness_diff=np.zeros([xxrange,36])
brightness_diff2=np.zeros([xxrange,36])


# for i in range(startx,xxrange+2):
    
startx=105  
# startx=120
endx=255    
# endx=125    


total_aur_bright=np.zeros([endx-startx])
total_me_bright=np.zeros([endx-startx])
total_po_bright=np.zeros([endx-startx])
total_eq_bright=np.zeros([endx-startx])

ii=np.arange(startx,endx)
ccml=(360-ii)/360*24  


for i in range(startx,endx):
# for i in range(startx+2,startx+3):

    

    cml=360-(90+i)    
    # cml=360-(90+i*20+10)   
    
    
    subsolarlat=0
    subsolarlongitude=cml#240.3#+180
    
    
    s_ssl=str(int(subsolarlongitude/10+1)*10)
    
    # print(s_ssl)
    
    
    v = open(planetary+'fieldline_tracing/fieldline_tracing_mapping_contours_kk2009ext_jrm09int_north_sslong'+s_ssl+'.txt')
    north_plotting = np.zeros([28,36,4])
    
    header = v.readline()
    
    for z in range(28):
        for y in range(36):
            line = v.readline()
            north_plotting[z,y,0]=line[0:12]
            north_plotting[z,y,1]=line[14:20]
            null1=line[21:29]
            if any(chr.isdigit() for chr in null1):
                north_plotting[z,y,2]=null1
            else:
                north_plotting[z,y,2]=np.nan
            null2=line[30:38]
            if any(chr.isdigit() for chr in null2):
                north_plotting[z,y,3]=null2
            else:
                north_plotting[z,y,3]=np.nan


    # ax1 = plt.subplot(1, 7, 8-i,projection=ccrs.Orthographic(360-polelong, polelat))

    mapfull=ab

    # mapmap=np.sum(mapfull[90+i*20:90+(i+1)*20,:,:],axis=0)
    mapmap=np.sum(mapfull[i:i+1,:,:],axis=0)
    mapmap[170:,:]=0 # because some local timers dont have emission here, so none should
    
    mapLT_map=np.rot90(np.sum(mapfull[90:270,:,:],axis=2),k=3)

    
    
    
    mapmapT=np.sum(mapfull[startx:endx,:,:],axis=0)
    mapmap=mapmap/np.mean(np.nonzero(mapmapT))
    mapmap_diff=mapmap/mapmapT
    # mapint=np.sum(a[75+i*30:75+(i+1)*30,:,:],axis=0)
    # mapcount=np.sum(b[75+i*30:75+(i+1)*30,:,:],axis=0)

    # mapmap=mapint/mapcount
    # mapmap=np.nan_to_num(mapmap)


    # mapmapT=totalmap

    # countmap=np.sum(b[90+i*20:90+(i+1)*20,:,:],axis=0)
    # totalmap3=np.sum(a[90+i*20:90+(i+1)*20,:,:],axis=0)
    countmap=np.sum(b[i:i+1,:,:],axis=0)
    totalmap3=np.sum(a[i:i+1,:,:],axis=0)
    
    # print(mapfull.shape,mapmap.shape,mapfull[75+i*30:75+(i+1)*30,:,:].shape,lon2d.shape)
    
    np_lon2=north_plotting[1,:,3]
    np_lat2=north_plotting[1,:,2]
    np_loc2=north_plotting[1,:,1]
   
    
    np_lon13=north_plotting[13,:,3]
    np_lat13=north_plotting[13,:,2]
    np_loc13=north_plotting[13,:,1]
    
    
    fit_xmo,fit_ymo=fit_fieldline_fixed(np_lon2,np_lat2) 
    fit_xmo2,fit_ymo2=fit_fieldline_fixed(np_lon13,np_lat13) 
    fity,fitx = fit_fieldline_scattered(ofl_lon,ofl_lat,sigma=6)
    # polar aurora should be included - we just want to measure total auroral brightness - teh shaping is necessary to remove pixel flaring
    # no- they should be done separately!
    
    mapmap2,lon2d2,lat2d2=masksubarray(lon2d,lat2d,mapmap,fit_xmo,fit_ymo)
    
    nonzerovals=np.nonzero(mapmap2)

    mapmap3=mapmap2[nonzerovals]
    lon2d3=lon2d2[nonzerovals]
    lat2d3=lat2d2[nonzerovals]


    magScale=True
    if magScale:

        

        param0=1120.572303190426
        param1=-43.641738021841704#g.parameters[0]

        B1=Bmag(lon2d3, lat2d3,radial=False)
        
        scaleint=B1*param1+param0
        # scaleint=scaleint#/np.mean(scaleint)

        mapmap3=mapmap3/np.reshape(scaleint,mapmap3.shape)

        B1=Bmag(lon2d, lat2d,radial=False)
        
        scaleint_eq=B1*param1+param0
        # scaleint=scaleint#/np.mean(scaleint)

        mapmap_eq=mapmap/np.reshape(scaleint_eq,mapmap.shape)

    mapmap4me,lon2d4me,lat2d4me=masksubarray(lon2d3,lat2d3,mapmap3,fit_xmo2,fit_ymo2,invert=True)
    mapmap4po,lon2d4po,lat2d4po=masksubarray(lon2d3,lat2d3,mapmap3,fitx,fity,invert=False)


    # plt.figure()

    # ax1 = plt.subplot(121)

    # ax1.tricontourf(360-lon2d4me, lat2d4me, mapmap4me,levels=128,cmap='afmhot') 
    # ax1.tricontour(360-lon2d4me, lat2d4me, mapmap4me,levels=128,cmap='afmhot',linewidths=1) 
    # ax1.set_xlim([90,230])
    # ax1.set_ylim([50,90])
    # ax2 = plt.subplot(122)
    # ax2.tricontourf(360-lon2d4po, lat2d4po, mapmap4po,levels=128,cmap='afmhot') 
    # ax2.tricontour(360-lon2d4po, lat2d4po, mapmap4po,levels=128,cmap='afmhot',linewidths=1) 
    # ax2.set_xlim([90,230])
    # ax2.set_ylim([50,90])
    

    
    total_eq_bright[i-startx]=np.median(np.append(mapmap_eq[90-40:90-25,:],mapmap_eq[90+15:90+30,:]))
    total_aur_bright[i-startx]=np.median(mapmap3)
    total_me_bright[i-startx]=np.median(mapmap4me)
    total_po_bright[i-startx]=np.median(mapmap4po)


    nplon1=north_plotting[3,:,3]
    nplat1=north_plotting[3,:,2]
    nplon2=north_plotting[13,:,3]
    nplat2=north_plotting[13,:,2]
    nplon3=north_plotting[23,:,3]
    nplat3=north_plotting[23,:,2]
    
    # ax2.scatter(nplon1, nplat1,transform=ccrs.PlateCarree(),c='k',vmax=maxint,vmin=minint) 



# In[9]



a_values=np.empty([3])
sigma_1d=np.empty([3])


g_init = models.Polynomial1D(degree=2)
#g_init = models.Gaussian1D()+ models.Polynomial1D(degree=2)
fit_g = fitting.LevMarLSQFitter()
# g = fit_g(g_init, ccml, total_me_bright,maxiter=1000)
g = fit_g(g_init, ccml, total_eq_bright,maxiter=1000)
# g = fit_g(g_init, ccml, total_eq_bright/np.max(total_eq_bright)*np.max(total_me_bright),maxiter=1000)


lt_param0=g.parameters[0]
lt_param1=g.parameters[1]


double_me = total_me_bright
double_me = np.append(double_me,np.flip(total_me_bright))
double_lt= np.append(ccml,ccml)



# plt.plot(ccml, total_me_bright,c='k',ls='--',alpha=0.5)
# plt.plot(ccml, g(ccml),c='C3',ls='--',alpha=0.5)



plt.show()

# plt.plot(np.abs(12-ccml), total_me_bright,c='k',ls='--',alpha=0.5)
# plt.plot(np.abs(12-ccml), g(ccml),c='C3',ls='--',alpha=0.5)


# plt.show()

g2 = fit_g(g_init, double_lt, double_me,maxiter=1000)
# plt.show()

# plt.figure()

# plt.scatter(ccml, total_me_bright,c='k',marker='.')
# plt.scatter(ccml, np.flip(total_me_bright),c='b',marker='.')
# plt.plot(ccml, g2(ccml),c='C3',ls='--',alpha=0.5,linewidth=3)
# plt.show()


noonlt=np.arange(100)/1000+12-0.05

# plt.plot(noonlt, g(noonlt),c='C3',ls='--',alpha=0.5,linewidth=3)


print(g.parameters)
print(g2.parameters)



# ax.plot(ccml,total_aur_bright)
fit_bright=g(ccml)
fit_bright2=g(np.arange(120)/10+6)
meanint=np.mean(fit_bright2)


i=1.5
# print((75+i*30)/360*24/24*360)


# plt.show()

ltrange_r1=np.arange(100)/50+7
ltrange_r2=np.arange(100)/50+9
ltrange_r3=np.arange(100)/50+11
ltrange_r4=np.arange(100)/50+13
ltrange_r5=np.arange(100)/50+15

# plt.plot(ltrange_r1,g(ltrange_r1)/meanint)
region1_scale=np.mean(g(ltrange_r1)/meanint)

# plt.plot(ltrange_r2,g(ltrange_r2)/meanint)
region2_scale=np.mean(g(ltrange_r2)/meanint)

# plt.plot(ltrange_r3,g(ltrange_r3)/meanint)
region3_scale=np.mean(g(ltrange_r3)/meanint)

# plt.plot(ltrange_r4,g(ltrange_r4)/meanint)
region4_scale=np.mean(g(ltrange_r4)/meanint)

# plt.plot(ltrange_r5,g(ltrange_r5)/meanint)
region5_scale=np.mean(g(ltrange_r5)/meanint)


# print(region1_scale,region2_scale,region3_scale,region4_scale,region5_scale)



# plt.scatter([8,10,12,14,16],[0.8828533281700559, 1.2033086368036778, 1.310723368138646, 1.2050975221749605, 0.8864310989126216])


LT_range=(np.arange(180)+90)/360*24
LT_scale=g(LT_range)/meanint

mapLT_map2=mapLT_map+0

for z in range(180): mapLT_map2[:,z]=mapLT_map2[:,z]/LT_scale[z]

# In[9]




tao10_lt= np.array([0.6927380920060449 ,1.4976335228606157 ,2.741562825090406  ,3.326941320257367  ,4.424525998695417  ,5.229421429549987  ,5.863581465980861  ,6.217247640144233  ,6.363592263935972  ,6.509936887727713  ,6.6562815115194525 ,6.802626135311193  ,6.948970759102933  ,7.095315382894674  ,7.278246162634348  ,7.534349254269893  ,7.863624657801308  ,8.266072373228594  ,8.668520088655878  ,9.144140115979035  ,9.847407335867118  ,10.671612126805321 ,11.522240252594809 ,12.41656850909989  ,13.132031114303953 ,13.717409609470913 ,14.193029636794066 ,14.632063508169285 ,15.03451122359657  ,15.363786627127986 ,15.693062030659402 ,16.058923590138754 ,16.38819899367017  ,16.680888241253648 ,16.973577488837126 ,17.266266736420608 ,17.55895598400409  ,17.851645231587568 ,18.290679102962788 ,19.002076579728193 ,19.818151113789074 ,20.660648982701247 ,21.501114287394085 ,22.251130484326755 ,22.97370706429847  ,23.668844027309238 ])
tao10_P= np.array([0.04975657151144641,     0.048050489576108135,    0.04465422581666478,    0.04462878563869144,    0.04123888192374159,    0.04120390167902821,    0.049253332990910836,    0.10800503900369551,    0.1614739330592,    0.2274760897943915,     0.28512273807645827,    0.33942718297727537,    0.38704722111559253,    0.42798285249141005,    0.47170206334047915,    0.5161143865383953,    0.5601336059862523,    0.6000554462698314,    0.6396430662152857,    0.678141290039587,    0.708701171329156,     0.7240302035673735,     0.7193977061595374,    0.6959634155522179,     0.6632530000513628,     0.6266755260689267,     0.5914423560147265,     0.5538234878421855,    0.5130311164685805,     0.4756676835830032,     0.4368002591758635,     0.3940041557846319,    0.3519616381653048,     0.31267802834663216,     0.26879888887874104,     0.2257553002561623,    0.18187616078827118,     0.14384587723756725,     0.10455590737440135,     0.07959772360629236,     0.06814306263890979,     0.06244327026453511,     0.058461087406143086,     0.05592183964217734,    0.0538015598092102,     0.05042914621661698])


tao10_P = tao10_P / np.max(tao10_P) * 1.35


fig, (z1,z2,z3) =plt.subplots(1,3,figsize=(12,5.5), gridspec_kw={'width_ratios': [7.5,10,10]})
plt.subplots_adjust(left=0.01,
                    bottom=0.1,
                    right=0.9,
                    top=0.9,
                    wspace=0.25,
                    hspace=0.15)
# mapmap4po,lon2d4po,lat2d4po

mapmapT2=mapmapT

# mapmapT2[170:,:]=0 # because some local timers dont have emission here, so none should

mapmap2,lon2d2,lat2d2=masksubarray(lon2d,lat2d,mapmapT2,fit_xmo,fit_ymo)

nonzerovals=np.nonzero(mapmap2)

mapmap3=mapmap2[nonzerovals]
lon2d3=lon2d2[nonzerovals]
lat2d3=lat2d2[nonzerovals]


magScale=True
if magScale:

    

    param0=1120.572303190426
    param1=-43.641738021841704#g.parameters[0]

    B1=Bmag(lon2d3, lat2d3,radial=False)
    
    scaleint=B1*param1+param0
    scaleint=scaleint#/np.mean(scaleint)

    mapmap3=mapmap3/np.reshape(scaleint,mapmap3.shape)

param0=1120.572303190426
param1=-43.641738021841704#g.parameters[0]

mapmap4me,lon2d4me,lat2d4me=masksubarray(lon2d3,lat2d3,mapmap3,fit_xmo2,fit_ymo2,invert=True)
mapmap4po,lon2d4po,lat2d4po=masksubarray(lon2d3,lat2d3,mapmap3,fitx,fity,invert=False)
B2=Bmag(lon2d, lat2d,radial=False)
scaleint2=B2*param1+param0

# scaleint2=np.reshape(scaleint2,mapmapT.shape)

subsolarlongitude=240
ax1 = plt.subplot(231,projection=ccrs.Orthographic(180, 80))
gl = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='black', alpha=0.5, linestyle='dotted', draw_labels=True)    
ax1.set_extent([150, 210, 52, 90], crs=ccrs.PlateCarree())
# cc=ax1.scatter(lon2d4po,lat2d4po,transform=ccrs.PlateCarree(),c=mapmap4po,vmin=0.4,vmax=1.2,cmap='afmhot',marker='.')
# ax1.scatter(lon2d4me,lat2d4me,transform=ccrs.PlateCarree(),c=mapmap4me,vmin=0.4,vmax=1.2,cmap='afmhot',marker='.')

gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
gl.ylocator = mticker.FixedLocator(np.arange(18)*10-90)




ax1.tricontourf(np.ravel(lon2d), np.ravel(lat2d), np.ravel(mapmapT/scaleint2),transform=ccrs.PlateCarree(),levels=128,cmap='afmhot',vmin=0.1) 
ax1.tricontour(np.ravel(lon2d), np.ravel(lat2d), np.ravel(mapmapT/scaleint2),transform=ccrs.PlateCarree(),levels=128,cmap='afmhot',vmin=0.1,linewidth=1) 



ax1.plot(north_plotting[1,:,3],north_plotting[1,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='k')
ax1.plot(north_plotting[1,:,3],north_plotting[1,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1.plot(north_plotting[13,:,3],north_plotting[13,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='k')
ax1.plot(north_plotting[13,:,3],north_plotting[13,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1.plot(np.arange(18)*10+90,np.zeros(18)+80,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='k')
ax1.plot(np.arange(18)*10+90,np.zeros(18)+80,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1.plot(fitx,fity,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='k')
ax1.plot(fitx,fity,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,c='grey',linestyle='--')

ax1.fill(np.append(np.roll(fitx,260),north_plotting[13,:,3]), np.append(np.roll(fity,260),north_plotting[13,:,2]),transform=ccrs.PlateCarree(), facecolor="none", hatch="///", edgecolor="k", linewidth=0.0,zorder=4)
ax1.fill(np.append(np.array([350,350,270,180,90,10,10]),north_plotting[1,:,3]), np.append(np.array([90,40,40,40,40,40,90]),north_plotting[1,:,2]),transform=ccrs.PlateCarree(), facecolor="none", hatch="//", edgecolor="k", linewidth=0.0,zorder=4)
ax1.fill(np.append(np.arange(70)*5+5,np.array([355,5])), np.append(np.zeros(70)+80,np.array([90,90])),transform=ccrs.PlateCarree(), facecolor="none", hatch='XX', edgecolor="grey", linewidth=0.0,zorder=5)

ax1.set_ylabel('Latitude')
ax1.set_xlabel('West Long.')


ax1.text(0.21, -0.08, '200$^{\circ}$W' ,transform=ax1.transAxes, ha="center", va="center")
ax1.text(0.5, -0.08, '180$^{\circ}$W' ,transform=ax1.transAxes, ha="center", va="center")
ax1.text(0.79, -0.08, '160$^{\circ}$W' ,transform=ax1.transAxes, ha="center", va="center")


gl.ylabels_right = False
gl.xlabels_top = False
gl.xlabels_bottom = False



scaled_equator_params = np.array([-1.01624303e-04  ,2.74302701e-05 ,-1.17873932e-06])



ax=plt.subplot(132)
ax.plot(ccml,total_me_bright/np.mean(total_me_bright)*np.mean(total_eq_bright)/meanint,c='C2')
ax.plot(ccml,total_po_bright/np.mean(total_me_bright)*np.mean(total_eq_bright)/meanint,c='C0')
ax.plot(ccml,total_eq_bright/meanint,alpha=0.5,c='C3')
ax.plot(ccml,total_eq_bright/meanint/np.mean(total_me_bright)*np.mean(total_eq_bright),c='C3')
ax.plot(tao10_lt,tao10_P,c='darkviolet',linestyle='dashdot')


tao_ccml=np.interp(ccml,tao10_lt,tao10_P)

# ax.plot(ccml,(ccml**2*scaled_equator_params[2]+ccml*scaled_equator_params[1]+scaled_equator_params[0])/meanint/np.mean(total_me_bright)*np.mean(total_eq_bright))
ax.set_ylim((0.,1.4))
ax.set_xlim((6,18))
ax.set_ylabel('Mean auroral brightness')
ax.set_xlabel('Planetary local time')


# ax.plot(ccml,total_eq_bright/meanint*12)

ax.plot(np.arange(120)/10+6, g(np.arange(120)/10+6)/meanint/np.mean(total_me_bright)*np.mean(total_eq_bright),c='hotpink',ls='--')
ax.plot(np.arange(120)/10+6, g(np.arange(120)/10+6)/meanint,c='hotpink',ls='--')



ax2=plt.subplot(133)
ax2.plot(ccml,total_me_bright/np.mean(total_me_bright)*np.mean(total_eq_bright)/fit_bright,c='C2')
ax2.plot(ccml,total_me_bright/np.mean(total_me_bright)/tao_ccml,c="darkviolet")
ax2.plot(ccml,total_po_bright/np.mean(total_me_bright)*np.mean(total_eq_bright)/fit_bright,c='C0')
ax2.plot(ccml,total_eq_bright/fit_bright,c='C3',alpha=0.5)
# ax2.plot(ccml,total_po_bright/np.mean(total_me_bright)*np.mean(total_eq_bright)/fit_bright,c='darkviolet')
ax2.set_ylim((0.5,1.5))
ax2.set_xlim((6,18))
ax2.set_ylabel('Ratio of local time enhancement')
ax2.set_xlabel('Planetary local time')



for i in range(1,7): ax.plot(np.array([(75+i*30)/360*24,(75+i*30)/360*24]),np.array([1.4,0.0]) ,c='k',ls='dotted',linewidth=0.5)
for i in range(1,7): ax2.plot(np.array([(75+i*30)/360*24,(75+i*30)/360*24]),np.array([1.5,0.5]) ,c='k',ls='dotted',linewidth=0.5)


mapLT_map_B=mapLT_map+0
mapLT_map_B[40:140,:]=mapLT_map_B[40:140,:]*5


ax1a=plt.subplot(234)
ax1a.imshow(mapLT_map_B, extent=[6,18,-90,90],cmap='afmhot',origin='lower')
ax1a.set_aspect(24/360)
ax1a.autoscale(False)
ax1a.fill(np.append(np.arange(70)*5+5,np.array([355,5]))/360*24, np.append(np.zeros(70)+80,np.array([90,90])), facecolor="none", hatch='XX', edgecolor="grey", linewidth=0.0)
ax1a.plot([6,18],[-40,-40],linewidth=1,alpha=0.5,c='k')
ax1a.plot([6,18],[-40,-40],linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1a.plot([6,18],[-25,-25],linewidth=1,alpha=0.5,c='k')
ax1a.plot([6,18],[-25,-25],linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1a.plot([6,18],[15,15],linewidth=1,alpha=0.5,c='k')
ax1a.plot([6,18],[15,15],linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1a.plot([6,18],[30,30],linewidth=1,alpha=0.5,c='k')
ax1a.plot([6,18],[30,30],linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1a.plot([6,18],[80,80],linewidth=1,alpha=0.5,c='k')
ax1a.plot([6,18],[80,80],linewidth=1,alpha=0.5,c='grey',linestyle='--')
ax1a.fill([6,18,18,6], [-25,-25,15,15], facecolor="none", hatch='XX', edgecolor="grey", linewidth=0.0)
ax1a.fill([6,7,7,6], [-90,-90,90,90], facecolor="none", hatch='XX', edgecolor="grey", linewidth=0.0)
ax1a.fill([18,17,17,18], [-90,-90,90,90], facecolor="none", hatch='XX', edgecolor="grey", linewidth=0.0)
ax1a.fill([6,18,18,6], [-40,-40,-90,-90], facecolor="none", hatch='XX', edgecolor="grey", linewidth=0.0)
ax1a.fill([6,18,18,6], [30,30,55,55], facecolor="none", hatch='XX', edgecolor="grey", linewidth=0.0)

# ax1a.plot([6,18],[-15,-15])
# ax1a.plot([6,18],[15,15])
# ax1a.plot([6,18],[30,30])
ax1a.set_xlim((6,18))

for i in range(1,7): ax1a.plot(np.array([(75+i*30)/360*24,(75+i*30)/360*24]),np.array([90,-90]) ,c='k',ls='dotted',linewidth=0.5)



bbox_props1= dict(boxstyle="circle,pad=0.15", fc="aliceblue", ec="C2", lw=1)
bbox_props2= dict(boxstyle="circle,pad=0.15", fc="linen", ec="C0", lw=1)
bbox_props3= dict(boxstyle="circle,pad=0.15", fc="mintcream", ec="C3", lw=1)
bbox_props3b= dict(boxstyle="circle,pad=0.15", fc="mintcream", ec="lightpink", lw=1)
bbox_props4= dict(boxstyle="circle,pad=0.15", fc="aliceblue", ec="black", lw=1)

ax1.text(180,70,'1', ha="center", va="center",bbox=bbox_props2, weight='bold',transform=ccrs.PlateCarree(),c='C0')
ax1.text(360-152,70,'2', ha="center", va="center",bbox=bbox_props1, weight='bold',transform=ccrs.PlateCarree(),c='C2')

ax1a.text(10,68,'1', ha="center", va="center",bbox=bbox_props2, weight='bold',c='C0')
ax1a.text(12,68,'+', ha="center", va="center", weight='bold',c='white')
ax1a.text(14,68,'2', ha="center", va="center",bbox=bbox_props1, weight='bold',c='C2')
ax1a.text(12.,21,'3', ha="center", va="center",bbox=bbox_props3, weight='bold',c='C3')
ax1a.text(12.,-33,'3', ha="center", va="center",bbox=bbox_props3, weight='bold',c='C3')

ax1a.set_ylabel('Latitude')
ax1a.set_xlabel('Planetary local time')

ax1a.set_xticks([6,9,12,15,18])

ax.text(17.5,0.86,'2', ha="center", va="center",bbox=bbox_props1, weight='bold',c='C2')
ax.text(17.5,0.95,'1', ha="center", va="center",bbox=bbox_props2, weight='bold',c='C0')
ax.text(17.5,0.07,'3', ha="center", va="center",bbox=bbox_props3, weight='bold',c='C3')
ax.text(17.5,0.073,'3', ha="center", va="center",bbox=bbox_props3b, weight='bold',c='C3',alpha=0.5)




ax1.text(360-260,75,'a', ha="center", va="center",bbox=bbox_props4, weight='bold',c='black',transform=ccrs.PlateCarree(),zorder=6)
ax1a.text(6.8,78,'b', ha="center", va="center",bbox=bbox_props4, weight='bold',c='black')
ax.text(6.5,1.35,'c', ha="center", va="center",bbox=bbox_props4, weight='bold',c='black')
ax2.text(6.5,1.47,'d', ha="center", va="center",bbox=bbox_props4, weight='bold',c='black')





fig.savefig('asd_lt.pdf', dpi=300, bbox_inches='tight', facecolor='white') 

# scaled_equator_params = [-1.01624303e-04  2.74302701e-05 -1.17873932e-06]

# ax1b=plt.subplot(4,3,10)
# ax1b.imshow(mapLT_map2, extent=[6,18,0,90],cmap='afmhot',origin='lower')
# ax1b.set_aspect(24/360)


# params = [-1.09858462e-04  2.79238756e-05 -1.16365027e-06] for 6-18 LT



# 0.07345374238775888


# In[7]


print(np.corrcoef(total_eq_bright,total_me_bright)) # 0.9357598
print(np.corrcoef(tao_ccml,total_me_bright)) # 0.80308569
print(np.corrcoef(total_eq_bright,total_po_bright)) # -0.17060778


plt.figure()
plt.plot(ccml,total_eq_bright)
# plt.plot(ccml,np.flip(total_eq_bright))
plt.plot(tao10_lt,tao10_P)
print(np.corrcoef(total_eq_bright,np.flip(total_eq_bright))) # 0.9357598



